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Abstract 

Michigan  Technological  University  (MTU),  under  contract  to  the  Air  Force  Office  of  Scientific  Re¬ 
search  (AFOSR),  has  been  investigating  strategies  for  overcoming  the  effects  of  scintillation  on  adaptive 
optical  imaging  and  laser  beam  projection  systems  which  must  operate  in  the  presence  of  extended  path 
turbulence.  The  approach  originally  proposed  under  this  contract  involved  using  two  deformable  mirrors 
to  fully  conjugate  the  received  or  transmitted  beam.  Our  investigations  into  this  concept  have  shown 
that  implementation  details  associated  with  the  use  of  real  adaptive  optical  hardware  will  likely  limit  the 
effectiveness  of  the  two  deformable  mirror  technique.  As  a  consequence,  an  alternative  approach  using 
one  deformable  mirror  was  developed,  and  its  performance  was  investigated.  This  technique  exploits  the 
fact  that,  under  strong  scintillation  conditions,  information  about  the  phase  of  the  beacon  field  is  incor¬ 
porated  in  the  intensity  of  the  beacon  field.  Hence,  deformable  mirror  commands  which  are  consistent 
with  all  of  the  intensity  measurements  available  to  an  adaptive  optical  system  are  processed  using  a  new 
nonlinear  optimization-based  algorithm.  Our  studies  of  the  performance  of  this  technique  have  shown 
that  under  conditions  of  saturated  scintillation  this  new  approach  significantly  outperforms  the  standard 
centroid-based  processing  of  the  Hartmann  sensor  measurements.  An  experiment  conducted  at  MTU  has 
demonstrated  the  ability  to  predict  a  laser  intensity  pattern  in  a  distant  target  plane  using  Hartmann 
sensor  measurements  of  the  field  reflected  from  a  deformable  mirror.  One  accepted  journal  article,  and 
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one  submitted  journal  article  have  resulted  from  this  work,  along  with  three  conference  papers.  In  ad¬ 
dition,  this  work  forms  the  basis  of  a  program  which  is  now  being  conducted  with  the  Airborne  Laser 
Technology  Division. 

1  Introduction. 

Atmospheric  turbulence  corrupts  propagating  optical  wave  fronts.  If  the  turbulence  lies  outside  the  near 
field  region  of  the  aperture  both  the  phase  and  the  amplitude  of  the  wave  are  corrupted  at  the  receiver 
[1,  2].  The  problem  of  correcting  for  the  phase  effects  of  atmospheric  turbulence  has  received  a  great  deal  of 
attention  in  the  literature  [3,  4].  Adaptive  optical  systems  which  can  effectively  conjugate  the  phase  of  a  wave 
front  corrupted  by  near  field  turbulence,  where  turbulence-induced  phase  effects  dominate  the  atmospheric 
degradation  of  performance,  are  becoming  mature  [5].  As  adaptive  optical  systems  continue  to  develop, 
interest  is  increasing  in  the  problem  area  of  projecting  laser  beams  and  imaging  over  longer,  and  possibly 
horizontal  paths,  where  both  phase  and  amplitude  effects  will  corrupt  the  propagating  wave  front  at  the 
receiver  or  target.  Amplitude  errors  in  a  wave  front  which  has  passed  through  a  turbulent  region  are  referred 
to  as  scintillation.  Scintillation  arises  in  optical  propagation  when  turbulence-induced  phase  aberrations 
in  one  region  are  converted  by  propagation  over  long  distances  to  amplitude  fluctuations  at  the  receiver 
or  target  for  the  system  [2].  Examples  of  systems  whose  performance  is  degraded  by  scintillation  include 
laser  weapon  systems  and  laser  communication  systems.  In  the  case  of  laser  weapons,  turbulence-induced 
scintillation  causes  amplitude  fluctuations  at  the  target  which  reduce  the  energy  density  that  can  be  applied 
at  the  desired  point.  Turbulence  induced  scintillation  degrades  the  performance  of  laser  communications 
systems  by  causing  signal  dropouts  which  reduce  the  signal-to-noise  ratio  [1]. 

The  amplitude  of  an  optical  wave  front  can  be  directly  controlled  using  absorbing  elements  such  as  a  liquid 
crystal  TV  [6].  Unfortunately,  absorbing  elements  are  inefficient  in  their  use  of  light,  and  are  inappropriate 
for  high  power  laser  applications  because  the  power  densities  involved  would  destroy  any  absorbing  element. 
Hence,  it  is  highly  desirable  to  accomplish  amplitude  control  using  entirely  reflective  elements  if  possible. 
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Reflective  deformable  mirrors  such  as  those  used  in  adaptive  optics  are  candidate  devices  for  amplitude 
control. 

An  approach  to  correcting  scintillation  effects  for  laser  beam  projection  and  imaging  systems  using  two 
deformable  mirrors  was  published  in  Ref.  [7],  which  is  an  extension  of  the  basic  concept  published  by  Kanev 
and  Lukin  in  Refs.  [8,  9,  10].  A  uniform  amplitude  beam  is  assumed  to  fall  on  the  first  deformable  mirror. 
The  location  of  the  second  deformable  mirror  is  chosen  to  lie  in  the  far  field  of  the  first  deformable  mirror, 
and  the  propagation  to  the  far  field  is  accomplished  by  use  of  a  Fourier  transforming  mirror.  The  problem  of 
determining  the  phase  adjustment  to  apply  to  the  first  deformable  mirror  is  not  linear,  and  cannot  be  solved 
in  closed  form  [11,  12].  However,  because  the  fields  falling  on  the  first  and  second  deformable  mirrors  form 
a  Fourier  transform  pair,  the  phase  of  the  first  deformable  mirror  can  be  obtained  using  an  approach  based 
on  the  Fienup  phase  retrieval  algorithm  [13].  The  second  deformable  mirror  is  used  to  correct  the  phase  of 
the  beam  to  achieve  field  conjugation  in  the  beam  transmitted  from  the  aperture.  This  mirror  configuration 
choice  enables  the  implementation  of  the  phase  retrieval  based  algorithm  for  determining  the  phase  of  the 
first  deformable  mirror  which  is  described  and  evaluated  here.  Because  the  phase  retrieval  algorithm  returns 
the  principal  value  phase  for  the  first  deformable  mirror,  phase  unwrapping  is  implemented  using  a  branch 
point  reconstruction  algorithm  [14]. 

The  results  published  in  Ref.  [7]  indicate  that  in  principle  the  two  deformable  mirror  scintillation  com¬ 
pensation  technique  can,  on  average,  improve  the  far  field  diflFraction  pattern  of  a  beam  passed  through 
extended  path  turbulence,  where  the  goal  is  to  maximize  the  energy  deposited  in  a  small  region  in  beacon 
plane.  However,  the  study  presented  in  Ref.  [7]  assumed  pointwise  control  of  the  deformable  mirrors.  This 
assumption  is  an  idealization  equivalent  to  assuming  that  the  deformable  mirrors  have  infinite  degrees  of 
freedom.  During  the  course  of  this  project  we  explored  the  effects  of  a  realistic  segmented  deformable  mirror 
model.  A  simulation  of  propagation  through  an  extended,  horizonatal,  constant  atmospheric  path  was 
used  to  evaluate  the  effectiveness  of  this  technique.  Two  compensation  techniques  were  implemented  to  cor¬ 
rect  an  outgoing  laser  beam:  (1)  conventional  phase-only  correction  based  on  branch  point  reconstruction  of 
the  beacon  field;  and  (2)  field  conjugation  implemented  by  the  two  deformable  mirror  technique  implemented 
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with  segmented  deformable  mirrors.  Since  the  focus  of  this  study  is  on  the  effects  of  realistic  deformable  mir¬ 
rors,  an  idealized  wave  front  sensor  is  assumed  which  can  measure  the  amplitude  of  the  incident  beacon  field 
and  local  phase  differences  on  the  spatial  sampling  scale  of  the  simulation.  The  results  of  our  study  indicate 
that  a  single  deformable  mirror  controlled  with  a  branch  point  reconstruction-based  approach  provides  equal 
or  superior  performance  compared  to  the  two  deformable  mirror  technique. 

As  a  consequence  of  this  finding  we  began  investigating  theoretical  avenues  to  process  all  of  the  intensity 
measurements  available  in  an  adaptive  optical  system  to  obtain  the  deformable  mirror  commands.  To 
understand  the  motivation  for  this  algorithm,  a  short  review  of  the  basic  theory  of  complex- valued  variables 
is  presented,  and  the  implications  adaptive  optical  systems  operating  in  the  presence  of  turbulence  are 
discussed.  Let  the  turbulence-corrupted  complex  optical  field  in  some  plane  be  represented  by  i7(x),  and 
adopt  the  standard  definition  of  the  principal  value  phase  ^(x)  as 

where  5R{(7(x)}  and  9{[/’(x)}  represent  the  real  and  imaginary  parts  of  the  complex  field,  respectively,  and 
X  represents  a  spatial  location  in  the  plane  of  interest.  Inspection  of  Eq.  (1)  shows  that  the  principal  value 
phase  is  defined  only  on  the  interval  (-7r,7r).  Hence,  in  an  optical  field  which  has  propagated  through  a 
turbulent  region  of  the  atmosphere  a  continuous  principal  value  phase  function  does  not  in  general  exist 
[15,  16].  There  are  two  sources  of  discontinuities  in  the  principal  value  phase:  (1)  so-called  “phase  wraps”, 
which  arise  from  the  phasor  representing  C/(x)  moving  from  the  second  to  the  third  quadrant  of  the  complex 
plane,  or  conversely;  and  (2)  so-called  “branch  points” ,  which  arise  from  the  presence  of  zeros  in  the  complex 
field.  When  plotted,  the  principal  value  phase  of  complex  fields  which  have  no  branch  points  may  appear 
discontinuous  due  to  phase  wraps,  but  in  fact  a  continuous  phase  map  with  values  that  lie  outside  the  range 
of  (— 7r,7r)  can  be  constructed  by  “unwrapping”  the  principal  value  phase  [14].  A  number  of  techniques  exist 
to  accomplish  this  unwrapping,  or  reconstruction  process.  These  techniques  are  based  on  computing  the 
discrete  gradient  of  the  principal  value  phase,  which  corresponds  to  computing  finite  phase  differences  on  a 


4 


sampled  grid  [3, 16],  and  the  applying  either  recursive  reconstruction,  such  as  used  in  reconstructing  the  phase 
of  the  Fourier  transform  of  the  image  in  the  Knox-Thompson  technique  of  speckle  imaging  method  [3],  least 
squares,  minimum  variance,  or  maximum  likelihood  techniques  for  fitting  smooth  elementary  functions  to 
the  phase  differences  [3] ,  or  by  computing  the  discrete  Laplacian  of  the  principal  value  phase  and  formulating 
a  solution  based  on  a  discrete  Poisson’s  equation  solution  [14], 

The  least  squares  reconstruction  paradigm  is  the  most  widely  used  in  adaptive  optical  systems  [3,  16]. 
However,  implicit  in  the  derivation  of  the  least  squares  reconstructor  is  the  assumption  that  there  exists  a 
continuous  phase.  As  discussed  in  Refs.  [15], [14],  and  [16],  when  there  are  branch  points  in  the  complex  field, 
a  continuous  phase  does  not  exist.  The  presence  of  branch  points  is  necessarily  associated  with  amplitude 
fluctuations,  or  scintillation,  in  the  optical  field.  Scintillation  strength  is  typically  characterized  by  the 
variance  of  the  logarithm  of  the  ratio  of  the  actual  amplitude  at  a  point  to  the  amplitude  of  the  field 
resulting  from  free  space  propagation  from  the  same  source,  denoted  by  .  Theoretical  expressions  for 
have  been  developed  using  a  weak  fluctuation  theory  [17,  18],  and  theoretically  obtained  values  for  are 
often  referred  to  as  the  “Rytov  variance” ,  which  we  will  denote  by  to  distinguish  the  theoretical  value 
for  the  variance  of  the  log  amplitude  fluctuations  from  the  experimentally  observed  value.  Under  conditions 
where  scintillation  is  negligible  is  small,  typically  on  the  order  of  (7^  ^  <  0.05,  and  and  ^  are 
typically  in  close  agreement.  As  the  path  length  L  or  the  strength  of  the  turbulence  (7^  increases,  ^ 
increases  steadily.  However,  it  has  been  observed  experimentally  [17]  and  in  simulations  [19,  20]  that  while 
cr^  ^  increases  steadily  as  or  L  increase,  the  observed  fluctuations  in  the  log  amplitude  saturate  at 
a  level  of  cr^  0.3  [17].  Under  conditions  of  ^  >  0.05  branch  points  generally  appear  in  the  turbulence 
corrupted  optical  field.  As  interest  grows  in  using  adaptive  optical  systems  to  image  and  propagate  laser 
beams  over  long  paths  or  through  strong  turbulence,  understanding  the  opportunities  for  improving  adaptive 
optical  system  performance  by  use  of  appropriate  wave  front  reconstructors  is  increasingly  important. 

A  new  class  of  wave  front  reconstruction  algorithms  have  been  developed  in  the  synthetic  aperture  radar 
signal  processing  community  which  can  sense  the  presence  of  branch  points  in  high  resolution  and  high 
signahto-noise  ratio  phase  difference  data,  and  provide  reconstructed  phase  maps  which  account  for  the 


5 


branch  points  and  associated  branch  cuts  [14].  One  branch  point  reconstruction  technique,  referred  to  as 
Goldstein’s  algorithm,  was  implemented  and  used  to  study  the  value  of  branch  point  reconstruction  in  laser 
beam  projection  systems.  Goldstein’s  algorithm  uses  a  two  step  process  to  reconstruct  a  wave  front:  (1)  it 
first  locates  branch  points  in  the  wave  front;  and  (2)  it  implements  a  path-dependent  phase  reconstruction 
similar  to  conventional  recursive  reconstructors  [3],  with  the  distinction  that  phase  discontinuities  arising 
from  branch  points  are  accounted  for  in  the  reconstruction.  When  branch  points  are  present  in  the  wave 
front,  Goldstein’s  algorithm  computes  a  discontinuous  phase  function  which  when  “rewrapped”  exactly 
matches  the  input  wrapped  pha^e  in  the  absence  of  noise.  Unfortunately,  a  significant  practical  barrier  to 
implementing  Goldstein’s  algorithm  is  the  need  to  explicitly  find  the  branch  points  in  realistic  wave  front 
sensor  data  which  is  spatially  averaged  and  corrupted  with  measurement  noise.  These  problems  may  prevent 
implementing  branch  point  reconstructors  with  real  wave  front  sensor  data,  and  provide  motivation  for  the 
study  of  a  nonlinear  approach  for  obtaining  deformable  mirror  commands. 

In  this  report  we  describe  an  algorithm  for  wave  front  reconstruction  and  deformable  mirror  control 
which  attempts  to  provide  performance  which  is  nearly  equivalent  to  branch  point  reconstructors  working 
with  pointwise  phase  difference  data,  but  does  not  require  that  the  branch  points  be  located.  Rather,  the 
algorithm  is  based  on  the  premise  that  information  about  the  phase  of  the  beacon  field  is  encoded  in  the 
intensity  of  the  beacon  field,  the  wave  front  sensor  intensity  image,  and  a  conventional  image  formed  with 
light  from  the  beacon.  The  algorithm  developed  here  is  an  extension  of  the  algorithm  presented  in  Ref.  [21], 
where  it  was  shown  that  Hartmann  sensor  and  conventional  images  could  be  jointly  processed  to  increase 
the  dynamic  range  of  a  Hartmann  sensor.  An  optimal  set  of  weights  for  wave  front  reconstruction  basis 
functions  or  deformable  mirror  influence  functions  can  be  obtained  by  finding  the  set  of  weights  which  is 
maximally  consistent  with  a  measured  Hartmann  sensor  image  and  a  measured  conventional  image  when 
the  pupil  intensity  is  available  to  the  algorithm.  Three  measurements  are  required  for  this  algorithm:  (1) 
the  Hartmann  sensor  image;  (2)  a  conventional  image  formed  with  light  from  the  beacon  field;  and  (3)  a 
beacon  intensity  map.  The  phase  of  the  incident  field  is  estimated  by  a  linear  combination  of  influence 
functions  in  a  nonlinear  optimization-based  algorithm.  For  comparison  purposes  we  have  also  implemented 
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idealized  least  squares  and  a  branch  point  reconstructor  based  on  Goldstein’s  algorithm,  both  of  which  use 
the  output  of  an  idealized  pointwise  phase  difference  sensor,  and  a  conventional  least  squares  reconstructor 
based  on  centroid  processing  of  Hartmann  sensor  output.  We  have  found  that  the  algorithm  developed  here 
provides  significantly  improved  performance  over  a  conventional  least  squares  Hartmann  sensor  processing 
algorithm  under  conditions  of  strong  and  saturated  scintillation,  and  with  performance  almost  equivalent  to 
Goldstein’s  algorithm  in  a  laser  beam  projection  paradigm  with  finite  degree-of-freedom,  phase-only  adaptive 
optics  implemented  in  the  transmitter.  One  benefit  of  the  algorithm  presented  here  is  that  the  need  for  point- 
by-point  phase  differences  required  by  Goldstein’s  algorithm  is  eliminated,  making  this  algorithm  attractive 
for  systems  which  must  operate  in  the  presence  of  scintillation. 

The  remainder  of  this  report  is  organized  as  follows.  The  theoretical  basis  and  supporting  analysis  are 
presented  in  Section  2.  The  simulation  used  to  evaluate  performance  is  discussed  in  Section  3.  Theoretical 
and  simulation  results  are  presented  in  Section  4,  and  experimental  results  are  presented  in  Section  5. 
Conclusions  are  drawn  in  Section  6. 

2  Technical  analysis. 

In  this  section  the  theoretical  analysis  of  the  scintillation  compensation  concept  is  presented.  This  section 
begins  with  the  basic  phase-retrieval  based  approach  to  the  laser  beam  projection  problem  which  assumes 
the  ability  to  control  an  optical  wave  front  with  a  very  high  degree  of  freedom  deformable  mirror.  Since 
deformable  mirrors  with  such  a  large  number  of  degrees  of  freedom  are  currently  an  idealization,  scintillation 
control  with  modal  mirrors  was  then  investigated.  Use  of  modal  mirrors  requires  a  significantly  different 
approach  to  determining  the  deformable  mirror  commands  since  such  mirrors  cannot  be  controlled  on  a 
“point-wise”  basis  as  was  assumed  for  the  early  work.  The  approach  for  determining  the  modal  mirror 
commands  is  based  on  nonlinear  optimization.  This  section  concludes  with  the  development  of  the  nonlinear 
optimization  theory  for  using  modal  mirrors.  The  scintillation  control  studies  implemented  for  imaging  follow 
directly  from  the  laser  beam  projection  work,  and  as  a  consequence,  a  less  detailed  theoretical  treatment  is 
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presented  for  the  imaging  work. 

The  functional  block  diagram  for  the  scintillation  correction  system  for  laser  beam  projection  is  shown  in 
Fig.  1.  It  is  assumed  that  a  beacon  wave  front  is  available  from  the  object  which  carries  the  amplitude  and 
phase  corruption  induced  by  propagation  from  the  object  to  the  telescope,  and  hence  carries  information 
about  the  beacon  field.  The  mechanism  by  which  this  beacon  is  provided  is  not  treated  here,  though  we  note 
that  obtaining  a  suitable  beacon  is  a  serious  issue  affecting  the  practical  implementation  of  this  technique. 
The  intensity  and  phase  of  the  incident  beacon  wave  front  is  measured  with  a  suitable  device,  such  as  a 
Hartmann  sensor,  which  is  sensitive  to  both  the  local  gradient  of  the  phase  and  the  intensity  of  the  incident 
wave  front  [3].  The  wave  front  sensor  outputs  are  processed  to  compute  the  amplitude  and  the  phase 
of  the  field  in  the  telescope  aperture.  Phase  retrieval  routines  process  the  beacon  amplitude  and  phase 
measurements  to  determine  the  phase  deformation  to  apply  to  the  deformable  mirrors  DM\  and  DM2  to 
conjugate  the  field  due  to  the  beacon.  It  is  assumed  that  a  collimated  laser  beam  is  incident  on  DMi,  The 
deformable  mirror  computer  contains  the  phase  retrieval  algorithm  which  computes  the  phase  distortion  to 
apply  to  DM\.  The  goal  is  to  apply  a  phase  distortion  to  DMi  which  provides  a  close  match  to  the  desired 
amplitude  distribution  when  the  light  falls  on  DM2*  DM2  lies  in  the  Fraunhofer  diffraction  region  of  DMi 
due  to  use  of  the  Fourier  transforming  mirror  shown  in  Fig.  1,  and  DM2  is  assumed  to  be  optically  conjugate 
to  the  pupil  plane  of  the  beam  projection  system.  The  purpose  of  DM2  is  to  correct  the  phase  distribution 
of  the  outgoing  beam  so  that  the  field  is  approximately  conjugated  to  the  incident  field  from  the  beacon  as 
it  leaves  the  aperture.  Two  approaches  for  controlling  DMi  are  now  discussed:  pointwise  control  of  DMi, 
and  modal  control  of  DMi. 

2.1  Pointwise  control  of  DMi. 

Let  the  complex- valued  wave  front  incident  on  the  aperture  from  the  beacon  Ub{xf)  be  represented  by 

Ub{xf)  =  Ab{xf)  exp  ,  (2) 
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Incoming  Scintillated  Light 
From  Beacon 


Figure  1:  Block  diagram  for  scintillation  correction  system. 


where  x/  is  a  position  vector  in  the  telescope  pupil,  Ab{xf)  is  the  amplitude  of  the  beacon  field,  is 

the  phase  of  the  beacon  field,  and  the  tilde  is  used  to  indicate  a  complex- valued  function.  Let  the  amplitude 
and  phase  of  the  beam  transmitted  by  the  aperture  Ut{xf)  be  represented  by 


Ut{xf)  =  B{xf)  exp  {jipti^f)}  . 


(3) 


where  B{xf)  represents  the  amplitude  of  the  field  reflected  from  DM2,  and  'iptixf)  represents  the  phase  of 
the  field  reflected  from  DM2-  In  the  laser  beam  projection  context,  optical  field  conjugation  means  that 
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B{xf)  and  are  adjusted  so  that  [22] 


B{S})  a  Ai{x!),  (4) 

and 

V’t(®/)  =  -Mxf),  (5) 

where  the  notation  zi  oc  Z2  is  used  to  represent  the  “proportional  to”  relational  concept.  Note  that  in  writing 
Eq.  (4)  we  have  assumed  that  the  amplitude  of  B{xf)  is  small  enough  so  that  there  are  no  nonlinear  effects 
such  as  blooming  [8,  23] 

The  amplitude  of  the  field  illuminating  DM2  is  established  by  applying  a  phase  deformation  to  DMi,  and 
projecting  the  resulting  field  into  the  Fraunhofer  diffraction  region  using  the  Fourier  transforming  mirror 
[24].  The  field  falling  on  DMi  is  assumed  to  be  a  well-collimated  laser  beam  of  amplitude  A.  The  field 
leaving  DMi  Ui{x)  has  amplitude  A  and  a  phase  profile  imparted  by  DMi 

Ui  (x)  =  A  exp  {jipi  (x) }  .  (6) 

Letting  the  field  illuminating  DM2  be  represented  by  1^2 (^/) 5  we  obtain 

U^ixf)  =  FT  [Ui (f)]  =  B{xf) exp  ,  (7) 

where  ^  is  a  spatial  coordinate  in  the  plane  of  DMi,  FT[‘]  represents  the  Fourier  transformation  operator, 
Xf  —  fXfi  where  /  is  the  independent  spatial  frequency  variable  which  arises  from  the  Fourier  transform 
operation  of  Eq.  (7)  and  fi  is  the  focal  length  of  the  Fourier  transforming  mirror,  'ipi  {x)  represents  the  phase 
distortion  applied  to  DMi,  and  represents  the  phgtse  of  the  wave  front  falling  on  DM2.  The  purpose 

of  DM2  is  to  apply  a  phase  distortion  to  U2{xf)  such  that  the  complex  transmitted  field  has  phase  equal  to 


The  relationship  between  the  phase  'ipi{x)  applied  to  DMi  and  the  amplitude  B{xf)  which  is  incident 
on  DM2  is  not  linear,  and  no  squared  error  metric  can  be  defined  which  would  yield  a  term  linear  in 
when  differentiated  with  respect  to  *01  Hence,  an  iterative  approach  to  determining  the  optimal  0i(f) 
is  adopted  here  which  is  based  on  the  phase  retrieval  algorithm  [12,  13].  A  functional  block  diagram  of  the 
algorithm  used  here  is  shown  in  Fig.  2. 

The  starting  point  for  the  algorithm  is  shown  in  the  upper  left-hand  corner  of  Fig.  2.  A  plane  wave  of 
amplitude  A,  which  physically  corresponds  to  the  amplitude  of  a  well-collimated  laser  beam  falling  on  DMi, 
is  associated  with  a  random  guess  for  0i(x),  the  phase  of  DMi.  The  field  leaving  J9Mi,  Aexp{j0i(f)}, 
is  then  projected  into  its  Fraunhofer  diffraction  region  by  use  of  a  fast  Fourier  transform  (FFT),  which 
physically  corresponds  to  the  function  performed  by  the  Fourier  transforming  mirror  shown  in  Fig.  1,  The 
outcome  of  this  step  is  represented  by  exp{j0(x/)},  and  physically  corresponds  to  the  field  which 

would  fall  on  DM2  due  to  the  particular  0i(x).  The  Fourier  domain  constraint  that  the  field  falling  on  DM2 
should  satisfy  Eq.  (4)  is  enforced  by  substituting  the  desired  field  amplitude  B'{xf)  for  B{xf)^  where  B'{xf) 
is  given  by 

{B  (^/))norm  ^  ('^*>(^/))norm  *  (^) 

The  subscript  “norm”  in  Eq.  (8)  is  used  to  indicate  that  the  terms  on  both  sides  of  the  equality  symbol  have 
been  normalized  to  have  unit  energy.  The  quantity  B'{xf)  exp{j0(f/)}  is  then  projected  back  into  the  plane 
of  DMi  by  use  of  an  inverse  FFT  to  obtain  A'(x)  exp{j0i(f)},  where  the  0i(^)  term  is  taken  to  be  the 
“next”  step  in  determining  the  ultimate  figure  of  DMi .  The  uniform  field  amplitude  A  is  then  substituted 
for  A'(x),  returning  the  algorithm  to  the  upper  left-hand  corner  of  Fig.  2.  This  process  is  repeated  until  the 
stopping  criterion  is  satisfied. 

The  stopping  criterion  used  here  is  based  on  the  convergence  of  a  measure  of  the  squared  error 
between  the  normalized  amplitude  of  the  field  falling  on  DM2^  (-^(^/))norm’  normalized  beacon 
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starting  Point 


^exp[yv^(  X  )  ] 


FT[.] 


B{  x^)exp[M  x^)  ] 


Space  Domain  Constraint  Fourier  Domain  Constraint 

A\x)  a  constant  amplitude  B{  Xj)  B\  jy)  =  [A  ^  x^)] 

Figure  2:  Block  diagram  of  phase  retrieval  technique  used  for  controlling  the  amplitude  of  the  field  falling 
on  DM2. 


field  amplitude  (^6  (^/))  norm  given  by 

~  J  [K^*’(^/))norml  “  K^(^/))norm|]  (^) 

Letting  a{k)  represent  the  value  of  a  obtained  in  the  iteration,  the  stopping  criterion  is  that  the  k^^ 
-01  (x)  is  taken  as  the  final  figure  of  DMi  when  the  condition 


\a{k)  -  a{k  -  1)\ 

CT{k  -  1) 


(10) 


is  satisfied.  We  note  that  this  stopping  criterion  does  not  guarantee  that  a  global  minimum  error  has  been 
found.  Rather,  this  stopping  criterion  indicates  that  0i(x)  has  stopped  changing  from  iteration  to  iteration. 
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2.2  Modal  control  of  DMi  and  DM2. 


The  derivation  presented  in  Section  2.1  assumes  that  the  phase  of  a  wave  front  can  be  controlled  on  a  pointwise 
basis.  This  level  of  control  is  an  idealization  which  cannot  now  be  realized  in  practice.  The  phase  retrieval 
algorithm  provides  the  principal  value,  or  “wrapped”  phase  (x)  due  to  the  use  of  Fourier  transforms  in 
the  calculation.  The  wrapped  phase  contains  discontinuities  due  to  the  complex  phasor  representation  of 
exp{'0i(x)}  moving  from  the  second  to  the  third  quadrant  of  the  complex  plane,  and  because  the  full  field 
propagated  backwards  from  the  DM2  plane  to  the  DMi  plane  after  the  constraints  are  imposed  in  general 
contains  branch  points.  Because  of  these  discontinuities,  fitting  a  deformable  mirror  directly  to  the  'ipi{x) 
which  is  output  by  the  phase  retrieval  algorithm  will  have  large  fitting  errors.  To  reduce  the  fitting  error 
we  have  implemented  a  branch  point  reconstruction  algorithm  known  as  Goldstein’s  algorithm  [14]  which 
we  use  to  “unwrap”  We  denote  the  unwrapped  phase  by  The  unwrapped  phase  will 

still  contain  discontinuities,  but  we  have  found  that  effects  due  to  fitting  error  are  much  smaller  when  phase 
unwrapping  is  implemented. 

The  segmented  mirror  was  specified  by  the  center-to-center  spacing  of  the  square  actuators,  and  it  was 
assumed  that  there  was  no  dead  space  between  actuators.  In  the  simulation  there  were  never  fewer  than 
3x3  pixels  used  to  model  individual  actuators.  The  positions  of  the  individual  actuators  was  determined 
by  assigning  the  actuator  to  have  the  same  value  as  at  the  location  of  the  center  location  of  the 

actuator.  Since  the  segmented  mirror  actuators  do  not  spatially  overlap  the  individual  influence  functions 
are  orthogonal,  and  hence,  no  additional  processing  is  required  to  determine  the  figure  for  DMi, 

The  field  falling  on  DMi  was  modeled  as  being  a  uniform  amplitude  plane  wave,  such  as  might  be 
obtained  from  a  well-collimated  laser  beam  with  amplitude  A.  The  field  falling  on  DM2,  denoted  in  by 
U2{xf),  was  calculated  using  Eq.  (7),  using  an  FFT  to  accomplish  the  Fourier  transform  operation.  The 
phase  of  i^M2,  denoted  by  must  compensate  both  for  the  phase  of  the  field  arriving  from  DMi^ 

represented  by  (f>{o^f)  in  Eq.  (7),  and  for  the  phase  of  the  atmosphere  so  that  in  principle 

^2iXf)  +  H^f)  =  (11) 
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In  practice,  we  have  found  that  the  field  falling  on  DM2  also  contains  branch  points,  and  hence  it  is  necessary 
to  unwrap  the  phase  which  we  wish  to  apply  to  DM2.  The  unwrapped  phase  for  DM2  found 

by  providing  the  wrapped  phase  given  by  arg{exp{-“0(rc/)}exp{-'0fc(aJ/)}}  to  Goldstein’s  algorithm.  A 
segmented  deformable  mirror  was  also  fit  to  obtain  the  actual  figure  for  DM2  The 

outging  field  was  calculated  by  computing 

Ut{xf)  =  U2{xf)exp{'ip2{x"f)}  ^  (12) 

2,3  A  new  algorithm  for  controlling  a  single  deformable  mirror. 

We  consider  wave  front  reconstruction  in  the  context  of  a  horizontal  path  adaptive  optical  laser  beam 
projection  system.  The  basic  geometry  is  shown  in  Fig.  3.  We  have  modeled  propagation  through  a  three 
dimensional  volume  of  turbulence  by  using  a  layered  model  for  the  atmosphere  [3, 19].  A  point  source  beacon 
illuminates  the  turbulence  volume,  and  the  turbulence-corrupted  beacon  field  is  intercepted  by  the  aperture. 
An  ideal  lens  is  placed  in  the  aperture  of  the  system  and  has  focal  length  equal  to  the  path  length  so  that  an 
outgoing  plane  wave  would  be  focused  at  the  beacon  location  in  the  absence  of  turbulence.  We  shall  refer  to 
the  complex  field  emanating  from  the  beacon  that  is  intercepted  by  the  telescope  pupil  as  the  beacon  field, 
and  denote  it  by  Ub{^p)  =  ^b(xp)  exp{j^B(xp)},  with  intensity  /b(xp)  =  |I7b(xp)|^,  where  xp  is  a  two 
dimensional  coordinate  in  the  pupil  plane.  A  single  deformable  mirror  and  a  Hartmann  wave  front  sensor  are 
placed  optically  conjugate  to  the  receiver  aperture.  In  addition  to  the  Hartmann  sensor  measurement,  it  is 
assumed  that  the  intensity  of  the  beacon  field  in  the  receiver  plane  can  be  measured,  and  that  a  conventional 
image  is  also  formed  with  the  beacon  field.  The  outgoing  beam  initially  has  uniform  intensity,  and  it  is  focused 
so  that  in  the  absence  of  turbulence  an  Airy  disk  would  be  present  at  the  beacon  location.  The  phase  of  the 
outgoing  beam  is  modified  using  a  simulated  deformable  mirror  whose  figure  is  adjusted  by  four  techniques 
to  perform  compensation  on  the  outgoing  beam:  (1)  a  conventional  Hartmann  sensor  centroid-based  least 
squares  reconstruction  algorithm  [3];  (2)  a  branch  point  reconstructor  known  as  Goldstein’s  algorithm  [14] 
which  requires  as  input  idealized  (and  in  practice  unobtainable)  pointwise  phase  differences;  (3)  a  least 


14 


Phase  Control 
System 


Beacon 


Random 

Phase 

Screens 


A 


-4- 

4- 

4 


vyu 

/ 

Receive/Transmit 

Aperture 


Outgoing 

Beam 


Figure  3:  Geometry  for  laser  beam  projection  system. 


squares  reconstructor  which  also  requires  as  input  idealized  (and  in  practice  unobtainable)  pointwise  phase 
differences;  and  (4)  the  new  nonlinear  optimization-based  algorithm  described  here.  It  should  be  noted  that 
the  branch  point  and  idealized  least  squares  reconstructions  would  not  be  obtainable  in  practice,  but  are 
computed  to  serve  as  comparative  measures  of  performance.  For  example,  the  performance  of  the  idealized 
branch  point  reconstructor  serves  as  an  upper  bound  on  the  performance  of  the  nonlinear  optimization- 
based  technique  developed  here,  and  the  idealized  least  squares  reconstructor  serves  as  an  upper  bound 
on  the  performance  of  the  Hartmann  sensor-based  least  squares  reconstruction.  The  compensated  beam  is 
propagated  back  through  the  turbulent  volume  to  the  beacon  plane,  where  statistical  quantities  of  interest 
are  accumulated  to  evaluate  the  performance  of  the  various  compensation  approaches. 

Since  we  are  interested  in  studying  the  effects  of  scintillation  on  wave  front  sensing  in  extended  path 
turbulence  we  have  modeled  a  100  km  path  using  10  equally  spaced  turbulent  layers,  with  the  first  layer  10 
km  from  the  beacon,  and  the  final  layer  lying  in  the  aperture.  The  ten  layer  model  is  justified  in  this  case 
by  the  close  match  between  the  Rytov  variance  computed  for  spherical  wave  propagation  through  a  region 
of  constant  given  by  [17] 


<^l,R  =  0.56A:'^/®  ^  dzC^iz)  (^)  ^  {L  -  zf/^, 


(13) 
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and  the  Rytov  variance  calculated  by  assuming  the  turbulence  lies  in  layers  given  by 


n=l  ^  ' 


(14) 


where  A;;,  =  10  km  is  the  layer  thickness.  For  example,  when  (7^  =  1  x  10  =  0.158,  while 


2  =  0.155. 

xJ 


layer 


The  phase  of  an  arbitrary  random  screen  is  denoted  by  ^fi(x),  where  x  is  a  coordinate  in  the  random 
screen  plane,  so  that  the  field  immediately  after  a  random  screen  has  the  form 


Ur{x)  =  [/yi(x)exp{j^i?(x)}, 


(15) 


where  Ua{x.)  is  the  field  incident  on  the  random  screen.  For  each  screen  the  field  Ur{x)  lies  in  the  Fresnel 
region  of  the  next  screen  or,  in  the  case  of  the  last  screen,  in  the  receive/transmit  aperture  [24].  An  angular 
spectrum  propagator  was  used  to  perform  the  near  field  propagations  [7].  The  beacon  field  is  passed  through 
a  simulated  lens  which  would  be  perfectly  focused  on  the  beacon  in  the  absence  of  turbulence,  and  thus 
having  complex  transmittance  given  by 

ti{xp)  =  exp  |-^^(|xi’l^)|  (16) 

where  the  total  path  length  L  =  100  km,  and  xp  is  a  two  dimensional  coordinate  in  the  pupil  plane.  Phase 
only  corrections  of  an  outgoing,  initially  planar  beam  are  calculated  in  the  aperture  plane,  and  the  corrected 
beam  is  passed  through  the  lens  and  propagated  in  reverse  through  the  phase  screens  to  the  beacon  plane. 
In  the  beacon  plane  the  intensity  distribution  of  the  various  correction  schemes  is  calculated,  and  statistical 
quantities  needed  to  compute  the  mean  and  the  variance  of  the  intensity  distribution  are  accumulated. 
These  intensity  statistics  form  the  basis  for  comparing  the  least  squares  and  branch  point  correction  schemes 
investigated  here. 

We  have  implemented  a  principal  phase  difference-sensitive  wave  front  sensor  with  resolution  equal  to  the 
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sample  spacing  in  the  simulation  to  provide  idealized  wave  front  sensor  input  into  the  idealized  branch  point 
and  least  squares  reconstructors.  The  wave  front  sensor  measures  principal  value  phase  differences  (or  some 
quantity  proportional  to  the  principal  value  phase  differences)  on  a  sampled  grid  in  xp-space  denoted  by 
the  discrete  integer  coordinates  with  0  <  i  <  M  —  1  and  0  <  j  <  N  -1.  Horizontal  phase  differences 
and  vertical  phase  differences  are  assumed  to  be  known  modulo  27r  on  a  discrete  grid  with 

locations  designated  by  The  phase  differences  are  defined  by  [14] 

A^{i,j)  =  arg[exp{j  +  l,i)  -  ipBiiJ))}] ,  (17) 

and 

A^{i,j)  =  arg[exp{j  (V’b(*,  j  +  1)  -  Mhi))}] .  (18) 

where  arg[-]  is  an  operator  which  takes  the  phase  of  a  complex-valued  quantity.  As  can  be  seen  by  Eqs.  (17) 
and  (18),  the  phase  differences  obtained  in  this  fashion  are  the  principal  value  phase  differences  which  are 
the  typical  output  of  wave  front  sensors  used  in  adaptive  optical  systems  [16].  Branch  points  in  the  field 
i7p(x)  will  rarely  (if  ever)  fall  directly  on  a  sample  point  [16,  14],  and  hence,  searching  for  a  zero  in  the 
incident  field  is  a  poor  strategy  for  locating  branch  points.  Rather,  the  phase  differences  are  summed  around 
all  possible  2x2  collections  of  samples  using 

S(i,i)  =  -A^{iJ)  -  Ay{i  +  1,  j)  +  A"(z,i  +  1)  -h  Ay{iJ).  (19) 

When  =  0  no  branch  point  is  enclosed,  and  no  phase  discontinuity  exists  at  (i,i).  However,  when 

=  ±27r,  a  branch  point  is  enclosed  by  the  loop,  and  it  is  expected  that  a  discontinuity  in  the 
reconstructed  phase  will  originate  at  this  point.  The  location  and  “polarity”,  or  sign  of  the  2ir  term  are 
noted  and  used  in  branch  point  reconstruction  algorithms.  The  branch  point  reconstruction  algorithm  used 
in  this  paper  is  generally  referred  to  as  the  Goldstein  algorithm  (though  the  seminal  paper  is  by  Goldstein, 
Zebker  and  Werner),  which  has  been  fully  discussed  elsewhere  [14].  Here  we  confine  ourselves  to  qualitatively 
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summarizing  the  algorithm  and  providing  an  example  of  its  performance.  Wrapped  phase  differences  are 
input  to  Goldstein’s  algorithm.  Branch  points  are  located  as  described  in  the  discussion  of  Eq.  (19),  and 
their  polarity  is  noted.  A  heuristic  approach  based  on  the  region  growing  concept  of  image  processing  [25] 
is  used  to  connect  the  nearest  pairs  of  positive  and  negative  branch  points.  Branch  cuts  are  constructed 
along  the  paths  defined  by  pairs  of  positive  and  negative  branch  points.  A  phase  reconstruction  approach 
similar  to  the  recursive  algorithm  is  then  implemented,  but  with  the  significant  exception  that  no  path  is 
allowed  to  encircle  a  branch  point  or  cross  a  branch  cut  As  a  direct  consequence,  the  reconstructed  phase  is 
discontinuous  at  the  branch  cuts.  The  unwrapped  phase  reconstructed  with  the  branch  point  reconstructor 
is  denoted  by  0gp(z,  j).  The  Goldstein  algorithm  is  known  to  provide  accurate  results  except  in  regions 
which  are  completely  enclosed  by  branch  cuts.  However,  in  our  experience,  the  Goldstein  algorithm  has  not 
failed  when  applied  to  phase  functions  created  by  simulating  laser  propagation  through  the  atmosphere. 

We  now  present  an  example  of  the  performance  of  the  Goldstein  algorithm.  Figure  4(a)  shows  the 
wrapped  phase  arg[exp{j?/;B(i,  j)}]  of  an  incident  beacon  field  simulated  as  described  in  Section  3.  This 
function  was  input  to  the  Goldstein  algorithm,  and  the  unwrapped  output  0gp(z,  j)  is  shown  in  Fig.  4(b), 
where  the  detected  branch  points  are  shown.  Figure  4(b)  is  hard  to  interpret  since  the  dynamic  range  of  the 
unwrapped  phase  is  much  larger  than  the  dynamic  range  of  the  wrapped  phase.  Hence,  in  Fig.  4(c)  we  show 
the  rewrapped  phase  given  by  arg  [exp{j0pp(z,  j)}].  Detailed  comparison  of  Fig.  4(a)  and  Fig.  4(c)  shows 
that  the  Goldstein  algorithm  provides  an  exact  phase  unwrapping  in  this  case  to  within  an  insignificant 
piston  term. 

Real  deformable  mirrors  have  a  finite  number  of  degrees  of  freedom,  and  as  a  consequence  </>gp(i,  j)  can 
not  in  general  be  fit  perfectly  with  a  real  deformable  mirror.  Here  we  modeled  a  deformable  mirror  with 
influence  functions  ejfc(xp)  consisting  of  a  collection  of  two  dimensional  linear  spline  functions,  also  referred 
to  as  two  dimensional  triangle  functions,  placed  on  a  Cartesian  grid  with  user  selected  spacing  D,  represented 
by  [26] 

e*,(xp)  =  A  (^^^-  "-)  ,  (20) 
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Wrapped  lixidont  Phase 


(a) 

Unwrapped  Inddent  Phase,  +  «  Pos.  BP.  o  -  Neg  BP 


(b) 

Re-wrapped  BP  Phase 


(c) 


Figure  4:  Example  of  branch  point  reconstruction  by  the  Goldstein  algorithm  for  a  laser  beam  propagation 
problem:  (a)  the  wrapped  incident  phase  ipj;  (b)  the  unwrapped  phase  (^/output  by  the  Goldstein  algorithm 
with  positive  and  negative  branch  points  indicated  by  (4-)  and  (o),  respectively;  and  (c)  the  rewrapped  phase 
given  by  arg[expO>/(i,;)}]. 
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where  is  the  center  location  of  the  actuator  projected  into  pupil  space.  It  should  be  noted  that 
the  bases  of  linear  splines  defined  this  way  are  2D  wide,  so  that  with  the  arrangement  used  here  the  line 
segments  defining  where  one  influence  function  goes  to  zero  fall  directly  under  the  peaks  of  the  adjoining 
influcence  functions.  These  actuators  were  fit  to  03p(i,i)  using  a  three  step  process  [3]: 

1.  The  reconstructed  phase,  is  projected  onto  the  actuators  using 


(21) 


where  Pk  is  the  projection  of  (j)^p{i,j)  onto  the  actuator,  and  W{iJ)  is  a  pupil  function  defined 
to  be  unity  inside  the  pupil,  and  zero  outside  the  pupil.  This  calculation  is  of  necessity  performed  at 
the  resolution  of  the  simulation  since  the  values  of  (f)^p{i^j)  only  exist  on  the  sample  grid.  The  entries 
of  Pk  can  be  arranged  into  a  vector  we  denote  by  p. 

2.  Since  the  influence  functions  are  not  orthogonal  the  influence  function  interaction  matrix  R  must  be 
computed  with  elements  given  by  [3] 


Tki  =  j  d^pW{Kp)ek{xp)ei{-Kp). 


(22) 


In  our  implementation  Eq.  (22)  is  computed  with  high  accuracy  using  numerical  quadrature  to  provide 
for  a  stable  inverse  of  R.  The  weight  vector  c,  containing  the  values  of  Ck  for  the  influence  functions, 
is  then  computed  from 

c  =  R^^p-  (23) 

3.  Finally,  the  deformable  mirror  phase  fit  to  0Bp(^,  j)  is  reconstructed  using 


=  '^Ckek{i,j), 

k 


(24) 


where  the  superscript  DM  is  used  to  represent  the  phase  of  the  deformable  mirror  applied  to  the 
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outgoing  beam. 


The  basic  principle  of  the  idealized  least  squares  wave  front  reconstruction  is  straightforward.  We  assume 
that  the  pointwise  wave  front  sensor  has  as  its  output  a  vector  of  x  and  y- directed  phase  differences  denoted 
by  A(j).  Representing  the  vectorized  input  phase  on  a  sampled  grid  using  the  notation  0,  we  can  define  the 
relationship 

=  P0,  (25) 

where  P  is  a  matrix  containing  Ts,  -I’s,  and  O’s,  and  is  ordered  so  that  the  correct  combinations  of  the 
elements  of  (j)  are  used  to  compute  each  element  of  [16,  14].  The  least  squares  estimate  of  (j)  is  given  by 

^LMSE  =  {P'^P)''P"’^.  (26) 

and  the  elements  of  the  vector  ^livisE  rearranged  into  an  array  denoted  by  0lMSE(^>  •?  )•  deformable 

mirror  is  then  fit  to  ^^ing  the  process  described  in  Eqs.  (21)  to  (24)  to  obtain  the  deformable 

mirror  figure  associated  with  the  idealized  least  squares  reconstruction  principle,  least 

squares  phase  reconstruction  based  on  Eq.  (26)  is  quite  simple,  but  in  practice  sparse  matrix  techniques  may 
be  required  to  implement  Eq.  (26).  When  phase  differences  computed  from  a  continuous  input  phase  are 
used  in  this  algorithm  in  the  absence  of  noise  the  output  phase  </>LMSE(^’-?)  ^^ifches  the  input  phase  to 
within  an  insignificant  piston  term  [16]. 

Due  to  the  use  of  pointwise  phase  difference  data  the  deformable  mirror  phase  0gp  {i,j)  and 
are  idealizations  that  can  not  be  obtained  experimentally.  Rather,  Hartmann  sensor,  or  other  wave  front 
sensor  data  is  generally  used  to  obtain  deformable  mirror  commands  using  a  least  squares  fitting  criterion, 
A  Hartmann  sensor  consists  of  an  array  of  lenslets  placed  one  focal  length  in  front  of  a  camera,  as  shown  in 
Fig.  5  [3].  Light  falling  on  an  individual  lenslet  is  focused  to  a  spot  in  the  focal  plane  which  is  in  general 
is  displaced  from  the  optical  axis  of  the  lenslet  due  to  turbulence  effects.  The  centroid  of  the  intensity 
distribution  of  the  light  in  the  detector  plane  can  be  calculated  for  each  lenslet,  and  the  local  slope  of  the 
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wave  front  can  be  calculated  using 


s  =  (27) 

Ji 

where  Xc  is  the  two  dimensional  centroid  location,  and  fi  is  the  focal  length  of  the  lenslet.  The  least  squares 
criterion  minimizes  the  sum  square  of  the  error  between  the  measured  slopes,  and  the  slopes  which  would  arise 
from  a  linear  combination  of  actuator  influence  functions.  This  minimization  is  accomplished  by  defining 
the  error 

A  =  8 -He  (28) 


where  s  is  a  vector  containing  the  measured  slopes  from  all  of  the  subaperture,  c  is  a  vector  containing  all  of 
the  deformable  mirror  weights,  and  H  is  a,  Jacobian  matrix  with  the  entry  in  the  row  and  column 
given  by 


jj-  _ 


(29) 


Minimizing  the  squared  error  defined  by  A^A  yields  the  least  squares  set  of  actuator  commands  [3] 


^LSH  ^  ^ 


(30) 


The  least  squares  estimate  for  the  deformable  mirror  figure  is  obtained  by  substituting  the  values  of  into 
Eq.  (24)  to  obtain  deformable  mirror  figure  obtained  from  the  Hartmann  sensor  data  j)- 

The  new  algorithm  is  an  extension  of  the  algorithm  presented  in  Ref.  [21],  where  it  was  shown  that 
the  entire  Hartmann  sensor  image  and  a  conventional  image  could  be  jointly  processed  to  increase  the 
dynamic  range  of  the  Hartmann  sensor.  This  algorithm  uses  a  parameterized  form  for  the  deformable  mirror 
figure  given  by  Eq.  (24),  but  with  a  diflFerent  approach  to  calculating  the  weights  Ck-  We  seek  to  estimate 
c  from  measurements  of  the  wave  front  sensor  image  intensity  /vv(xw),  where  x^  is  a  two  dimensional 
coordinate  in  the  Hartmann  sensor  image  plane,  the  pupil  plane  intensity  /p(xp),  and  the  conventional 
image  plane  intensity  /c(xc))  where  xc  is  a  two  dimensional  coordinate  in  the  conventional  image  plane. 
The  optimization  is  accomplished  by  forming  an  objective  function  J(c),  and  minimizing  J(c)  by  the  choice 
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Figure  5:  Diagram  of  the  Hartmann  wave  front  sensor. 


of  c.  The  objective  function  used  here  is  given  by 


(ic(xc,c)^  ,  (31) 

^  -*  V/-*  ^ 


where  /^^(xwjc)  and  /c(xcjc)  represent  the  wave  front  sensor  and  conventional  images  computed  from  the 
current  guess  of  the  vector  c,  respectively.  The  objective  function  in  Eq.  (31)  is  an  extension  of  one  of  the 
objective  functions  used  in  phase  retrieval  [27].  Also,  the  availability  of  the  beacon  field  intensity  information 
/j5(xp)  provides  a  significant  extension  to  the  approach  used  in  [21]  due  to  the  fact  that  [/p(xp)]^^^  = 
\Ub{'x.p)\  can  be  used  in  computing  iW(xwjC)  and  JcCx/jc)  as  described  below.  The  objective  function 
in  Eq.  (31)  can  be  minimized  using  quasi-Newton  methods  [28,  29],  either  with  or  without  the  use  of  an 
analytic  gradient  vector  VJ,  with  element  given  by 


dJ{c) 

dck 


(32) 


If  no  analytic  gradient  is  used,  then  it  is  necessary  to  perturb  each  parameter  Ck  at  each  step  to  evaluate 
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both  VJ  and  the  Hessian  matrix  H  of  J(a)  with  element  given  by 


Hki 


5V(c) 

dckdci  ‘ 


(33) 


This  perturbation-based  approach  to  estimating  V  J  and  H  is  computationally  expensive.  Hence,  an  analytic 
gradient  V  J  was  calculated.  The  details  of  this  analysis  are  described  in  [21],  and  we  present  only  the  results 
of  the  analysis  here.  Note  that  the  objective  function  in  Eq.  (31)  has  two  terms  which  are  summed:  one 
due  solely  to  the  wave  front  sensor  image,  and  one  due  solely  to  the  conventional  image,  and  we  let  these 
two  contributions  to  the  objective  function  be  represented  by  Jw{c)  for  the  wave  front  sensor  channel,  and 
Jc{c)  for  the  conventional  image  channel.  The  gradients  of  these  two  terms  are  given  by 


dJwjc) 

dck 


=  2^  \^tHixp)ekiy.p)exp{j 


yv  ^  Iwj  Uwi^w,c)hwi'x.p,xw) 

^  \Iw{x.w,c)J 


(34) 


and 


dJcjc) 

dCk 


=  2Q 


l  xp  I  fc 


/  Jc(xc)  \ 
\Ic{x.i,c)J 


1/2 


^c(xc,c)/ic(xp,xc) 


(35) 


where  xp  represents  a  pupil  plane  coordinate,  xc  represents  a  conventional  image  plane  coordinate,  and 
xw  represents  a  Hartmann  sensor  image  plane  coordinate.  In  Eq.  (34)  the  term  t/f(xp)  is  the  complex 
transmittance  of  a  Hartmann  sensor  lenslet  array  with  square  lenslets  of  side  length  d  and  focal  length  fi 


tH(y-p)  =  Xexpjj^ixp  -Xc(5)p|rect  .  (36) 

where  =  2it/\,  X  is  the  mean  wavelength  of  the  light,  S  is  the  number  of  subapertures,  and  Xc(s)  is 
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the  center  of  the  subaperture,  the  term  hwi^p^'^w)  is  the  propagation  kernel  for  propagation  from  the 

back  of  the  lenslet  array  to  the  lenslet  focal  plane,  U^{xW)C)  is  the  conjugate  of  the  field  estimate  in  the 

«  2 

wave  front  sensor  detector  plane  associated  with  the  current  state  of  c,  and  /w(xiy,c)  =  Uwi^WiC)  .  The 

propagation  from  the  back  of  the  Hartmann  sensor  lenslet  array  to  its  detector  plane  is  implemented  with 

an  angular  spectrum  propagator  [21,  24],  and  the  field  propagated  is  given  by  |I7b(xp)|  exp{5]]jj.  CkekiiJ)}^ 

incorporating  information  about  the  amplitude  fluctuations  of  the  beacon  field  into  the  algorithm.  In  Eq.  (35) 

the  term  P{xp)  represents  the  physical  extent  of  the  pupil,  UQ{xcyC)  is  the  conjugate  of  the  conventional 

2 

image  field  estimate  associated  with  the  current  state  of  c,  Ici'x.i^c)  =  Uci^c^c)  ,  and  hc{y^p,^p)  is  the 
propagation  kernel  from  the  pupil  plane  to  the  conventional  image  plane,  and  the  amplitude  of  the  field 
propagated  is  given  by  |f7B(xp)|.  The  propagation  from  the  pupil  to  the  conventional  image  plane  can  be 
accomplished  computationally  using  Fourier  transform  techniques  [21,  24]. 

The  weights  obtained  from  the  optimization  process  are  used  in  Eq.  (24)  to  form  the  deformable  mirror 
figure  associated  with  the  nonlinear  optimization  algorithm  denoted  by  In  the  next  section  we 

describe  a  simulation  for  propagating  a  point  source  beacon  field  through  a  turbulent  region  to  a  receiving 
aperture,  applying  the  desired  correction  to  the  outgoing  beam,  and  propagating  the  fields  resulting  from 
the  three  phase  correction  paradigms  used  here  from  the  aperture  back  through  the  turbulence  and  on  to 
the  beacon  plane,  where  turbulence  compensation  results  can  be  compared. 

3  Simulation  description 

The  algorithm  developed  here  is  not  linear,  and  further,  we  wish  to  evaluate  its  performance  for  saturated 
scintillation  situations.  Hence,  a  simulation  was  used  to  evaluate  the  performance  of  the  new  algorithm.  The 
simulation  follows  the  basic  conceptual  form  shown  in  Fig.  3.  The  field  due  to  a  point  source  beacon  was 
calculated  in  the  plane  of  the  first  phase  screen  a  distance  10  km  away  from  the  beacon  [24],  The  beacon  field 
was  propagated  through  the  phase  screen  using  Eq.  (15).  At  each  step  the  field  which  was  passed  through 
the  phase  screen  was  propagated  10  km  to  the  next  phase  screen  using  an  angular  spectrum  propagator 
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[21,  24,  7],  and  the  process  was  repeated  until  the  final  screen  was  encountered  in  the  aperture  of  the  system. 
A  total  of  10  random  phase  screens  were  used  to  simulate  the  100  km  optical  path.  All  of  the  phase  screens 
and  propagations  were  conducted  in  256  x  256  arrays  with  1  cm  sample  spacing.  As  discussed  in  Ref.  [7], 
the  relationship  between  the  number  of  points  along  one  edge  the  fast  Fourier  transform  (FFT)  array  N, 
the  sample  spacing  Ax,  the  wavelength  A,  and  the  propagation  distance  z  is 


N> 


2Xz 
(Ax)^  * 


(37) 


For  the  case  used  here  with  A  =  987  nm,  z  =  10  km,  and  Ax  =  1  cm,  using  Eq.  (37)  yields  N  >  197.4,  so 
that  working  with  N  =  256  -  sized  arrays  is  satisfactory. 

Because  a  point  source  beacon  emanates  waves  in  all  directions  it  is  necessary  in  a  simulation  to  strongly 
attenuate  rays  which  are  well  outside  the  solid  angle  subtended  from  the  aperture  to  the  beacon  [19].  Atten¬ 
uating  these  rays  prevents  “wrap-around  error”  which  would  arise  in  an  FFT-based  propagation  simulation 
if  these  rays  were  allowed  to  propagate  unattenuated.  Our  strategy  for  attenuating  these  rays  is  derived 
from  Ref.  [19].  Over  a  region  which  is  1.5  times  the  diameter  of  the  circles  described  by  the  intersection  of 
the  solid  angle  subtended  from  the  aperture  to  the  beacon  and  the  planes  containing  the  phase  screens  the 
beam  is  allowed  to  propagate  unattenuated.  Outside  this  region  a  Gaussian  rolloff  is  used  to  attenuate  the 
field  with  e”^  width  of  25  samples.  For  the  propagation  from  the  aperture  to  the  beacon  plane  the  field  is 
converging,  and  as  a  result  it  is  not  necessary  to  attenuate  the  field  in  this  manner. 

The  phase  screens  were  generated  using  a  white  noise  filtering  approach  [19,  20,  30,  31].  Phase  screens 
possessing  the  von  Karman  spectrum  with  outer  scale  Lq  =  100  m  were  simulated  using  this  method  so  that 
ensemble  of  phase  screens  has  the  power  spectrum 

0m969eziCl 

where  /  is  a  spatial  frequency  with  units  of  inverse  meters,  and  zi  =  10  km  is  the  layer  thickness.  The 


26 


strength  of  each  screen  is  controlled  by  adjusting  The  phase  screens  were  created  in  2N  x  2N  arrays, 
and  the  center  N  x  N  section  was  extracted  to  avoid  problems  arising  from  the  fact  that  autocorrelation 
of  the  phase  is  an  even  function  [31].  While  this  approach  to  phase  screen  generation  has  been  widely  used 
[19,  20,  30,  31,  32],  this  approach  is  known  to  produce  phase  screens  which  are  affected  by  truncation  of 
the  power  spectrum  at  both  high  spatial  frequencies,  due  to  finite  sample  spacing  Ax,  and  at  low  spatial 
frequencies,  due  to  finite  array  size  N.  We  adopt  this  phase  screen  generation  technique  here  because  of  its 
speed  of  operation,  and  because  of  reasonable  agreement  between  theory  and  simulation  [19]. 

The  beacon  field  intercepted  by  the  aperture  was  passed  through  a  simulated  lens  with  complex  trans¬ 
mittance  given  by  Eq.  (16).  The  intensity  of  the  beacon  field  Ib{^p)  is  noted,  and  used  to  compute  the 
amplitude  of  the  beacon  field  |{7j3(xp)|,  which  is  used  as  the  field  amplitude  in  computed  propagations  to  the 
Hartmann  sensor  and  conventional  image  planes  conducted  in  the  nonlinear  optimization-based  algorithm. 
A  conventional  image  Ici^c)  and  a  Hartmann  sensor  image  Iw{^w)  are  computed  for  each  beacon  field, 
and  this  information  is  also  provided  to  the  nonlinear  optimization-based  algorithm  to  compute 
A  conventional  Hartmann  sensor  centroid-based  least  squares  wave  front  reconstructor  [3]  was  implemented 
on  the  wave  front  sensor  image  to  obtain  0lSH(^’-?)’  ^  simulated  outgoing  laser  beam  initially  carrying 

the  resulting  phase  was  propagated  to  the  beacon  plane  and  evaluated  to  provide  one  comparative  perfor¬ 
mance  measure  for  the  new  algorithm.  Point- wise  principal  value  phase  differences  were  also  computed  for 
the  beacon  field  using  Eqs.  (17)  and  (18),  which  are  input  to  the  branch  point  phase  reconstructor  to  obtain 
0Bp(^ji)  and  the  ideal  least  squares  reconstructor  to  obtain  0LMSE(^»-?)*  ^  deformable  mirror  was  fit  to 
<^gp(i,  j)  and  ^LMSE^^j-^)  obtain  («,  j)  and  these  phases  were  applied  to  the  out¬ 

going  beam  to  obtain  the  branch  point  and  ideal  least  squares  reconstruction  results.  It  should  be  noted  that 
these  pointwise  phase  differences  are  an  idealization  of  a  real  wave  front  sensor,  and  would  not  be  available 
in  practice.  Hence,  these  results  represent  an  upper  bound  on  performance  that  could  not  be  obtained.  We 
provide  these  results  to  obtain  another  comparative  measure  of  performance  for  the  new  algorithm. 

The  outcomes  of  the  four  phase  correction  schemes  were  associated  with  an  initially  planar  wave  in  the 
aperture  to  create  an  outgoing  field  which  was  propagated  to  the  beacon  plane  and  evaluated  for  comparison 
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to  the  performance  of  the  new  algorithm.  The  outgoing  fields  were  propagated  back  through  the  simulated 
lens  and  the  phase  screens,  in  reverse  order,  to  the  beacon  plane.  In  the  beacon  plane  statistical  quantities 
which  are  used  to  compute  the  mean  intensities  for  all  of  the  correction  approaches  studied.  In  the  absence 
of  turbulence,  the  intensity  pattern  in  the  beacon  plane  would  be  an  Airy  disk  [24],  but  with  turbulence 
and  the  various  correction  schemes  the  actual  beacon  plane  patterns  are  suboptimal.  Statistical  quantities 
required  to  compute  the  average  far  field  intensity  and  the  average  Strehl  ratio  were  accumulated,  where  the 
appropriate  definition  of  Strehl  ratio  SR  for  this  case  is 


SR  = 


Actual  on-axis  intensity 
On-axis  intensity  with  no  atmosphere* 


(39) 


We  also  computed  the  radial  average  of  the  intensity  patterns  in  the  beacon  plane.  We  now  detail  the 
physical  parameters  used  in  the  simulation. 

A  0.5  m  diameter  receive/transmit  aperture  was  modeled,  and  the  wavelength  of  operation  was  set  at  987 
nm.  Deformable  mirror  actuators  were  placed  on  a  Cartesian  grid  with  spacing  of  8  cm  in  the  plane  of  the 
receive/transmit  aperture,  providing  for  a  total  of  37  actuators.  Deformable  mirror  influence  functions  were 
implemented  as  described  in  Eq.  (20).  When  equal  weights  are  applied  to  all  actuators  in  this  configuration 
a  planar  deformable  mirror  figure  results.  The  wave  front  sensor  had  8  subapertures  across  its  diameter  for  a 
total  of  35  fully  illuminated  subapertures  within  the  pupil.  This  arrangement  of  actuators  and  subapertures 
provides  a  stable  solution  to  Eq.  (30).  When  projected  into  the  pupil  plane  the  subapertures  had  side  length 
equal  to  6  cm,  though  they  were  actually  modeled  as  having  physical  side  length  of  90  fim  and  focal  length 
of  3  mm.  The  sample  spacing  in  the  wave  front  sensor  lenslet  and  detector  planes  was  15  /xm,  satisfying  the 
sampling  criterion  for  the  wave  front  sensor  lenslet  array  [21] 

Axw<^,  (40) 


where  Axw  is  the  sample  spacing  in  the  wave  front  sensor  plane  since  Axw  =  15/im  <  21.9/im.  The  image 
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plane  for  each  subaperture  occupied  a  7  x  7  sample  region  behind  each  lenslet.  A  Monte  Carlo  simulation 
approach  was  used,  with  50  statistically  independent  propagations  used  to  obtain  the  results  presented. 

The  simulation  can  be  evaluated  by  comparing  the  theoretical  Rytov  variance  calculated  using 
Eq.  (13)  to  the  variance  of  the  log  amplitude  fluctuations  at  the  aperture  calculated  from  a  set  of  100 
simulated  propagations  through  the  turbulent  volume.  This  plot  is  shown  in  Fig.  6,  where  the  cr^  ^  and  the 
simulated  are  plotted  as  a  function  of  for  the  path.  Inspection  of  Fig.  6  shows  that  the  simulation 
consistently  provides  slightly  stronger  values  of  than  are  be  predicted  by  Eq.  (13)  when  <t^  <  0.3,  and 

the  cr^  provided  in  the  simulation  saturates,  as  expected  [17,  19],  at  «  0.35.  The  fact  that  simulated 
is  larger  than  for  low  values  of  most  likely  arises  from  the  attenuation  of  the  rays  implemented  to 
avoid  wrap-around  error  discussed  above,  and  the  limitations  of  the  phase  screen  generator,  also  discussed 
above. 


4  Theoretical  and  simulation  results. 

Results  of  this  study  are  presented  in  Figs.  7  and  8.  Figure  7  shows  radial  averages  of  the  far  field  intensity 
patterns  for  the  C^  =  3xl0“^’^  case,  where  the  peak  intensity  has  been  normalized  to  the  intensity  obtained 
by  the  branch  point  reconstruction  technique.  In  this  case  the  Rytov  variance  is  0.47,  and  saturated  log 
amplitude  fluctuations  were  obtained  in  the  simulation.  Inspection  of  Fig.  7  shows  that  the  branch  point 
reconstructor  performed  significantly  better  than  the  least  squares  reconstructor,  a  consistent  result  for  all 
cases  where  saturated  log  amplitude  fluctuations  arose. 

Figure  8  shows  SR  vs.  and  SR  vs.  the  Rytov  variance  for  all  of  the  simulation  runs  conducted. 
Inspection  of  Fig.  8  shows  that  for  low  values  of  Rytov  variance,  or  equivalently,  low  values  of  the  branch 
point  and  least  squares  recontructors  perform  comparably.  However,  in  the  region  where  the  log  amplitude 
fluctuations  are  saturated,  >  0.35,  branch  point  reconstruction  has  a  significant  advantage  in  Strehl  ratio 
over  least  squares  reconstruction.  These  results  demonstrate  the  importance  of  branch  point  reconstruction 
for  adaptive  optical  laser  beam  propagation  through  the  atmosphere  for  optical  paths  that  cause  saturated 
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Figure  6:  Theoretical  Rytov  variance,  and  simulated  cr^  as  a  function  of  (7^. 

log  amplitude  fluctuations.  Further,  these  results  show  that  fitting  errors  between  the  deformable  mirror 
and  the  discontinuous  branch  point  reconstructed  phase  do  not  dominate  the  errors  in  the  outgoing  beam. 

The  main  results  of  this  paper  are  presented  in  Figs.  9  and  10.  The  input  (7^  values,  the  associated 
Rytov  variances  and  the  average  number  of  objective  function  evaluations  required  convergence  Nj 

are  listed  in  Table  1.  In  the  legends  of  both  these  flgures  “BP”  refers  results  obtained  for  “LS” 

refers  to  results  obtained  for  j),  “WFS”  refers  to  results  obtained  with  and  “NL  OPT” 

refers  to  results  obtained  with  j).  Results  for  the  curves  labeled  BP  and  LS  are  unobtainable  in 

practice  because  they  require  that  pointwise  phase  differences  be  available,  and  these  curves  are  provided  for 
purposes  of  comparison.  The  curves  labeled  WFS  and  NL  OPT  are,  in  principle,  obtainable  in  practice  since 
they  can  be  computed  from  intensity  information  available  in  an  adaptive  optical  system.  The  branch  point 
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Radial  Average  Intensity  Patterns 


Figure  7:  Radial  averaged  far  field  intensity  patterns  for  the  case  (7^  =  3  x  10 

reconstruction  results  represent  the  upper  bound  on  performance  for  the  nonlinear  optimization  technique, 
and  the  idealized  least  squares  curves  represent  an  upper  bound  on  the  performance  of  the  Hartmann  sensor- 
based  least  squares  reconstruction  technique. 

Figure  9  shows  the  average  Strehl  ratio  results  for  the  four  correction  techniques  investigated.  For 
7  X  10“^®  <  <  1  X  10“^^  scintillation  is  weak,  and  using  conventional  Hartmann  sensor  processing 

provides  better  Strehl  ratio  performance  than  the  nonlinear  optimization  technique.  However,  as  cr^  ^ 
increases,  the  nonlinear  optimization  technique  begins  to  outperform  the  conventional  Hartmann  sensor 
processing  approach.  The  onset  of  saturated  scintillation  occurs  in  the  vicinity  of  =  2  x  10”^^  for  the 
geometry  studied  here.  In  the  regime  >  2  x  the  nonlinear  optimization  technique  provides  Strehl 
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Strehl  Ratio 


Strehl  Ratio  for  Various  Reconstruction  Techniques 


(a) 


Rytov  Variance 

(b) 


Figure  8:  Strehl  ratio  results  for  branch  point  and  least  squares  reconstruction,  and  for  no  compensation: 
(a)  Strehl  ratio  vs.  and  (b)  Strehl  ratio  vs.  Rytov  variance 


Table  1:  Listing  of  and  Nj  for  the  plots  shown  in  Figs.  9  and  10. 


c'i 

- n - 

Ni 

7  X  10“^* 

0.11 

44.5 

1  X  10-1'^ 

0.16 

49.9 

2  X  10-1’’ 

0.32 

53.6 

3  X  10-1^ 

0.47 

57.5 

4  X 

0.63 

56.8 

5  X  10-1^ 

0.79 

58.5 

6  X  10-1^ 

0.95 

53.6 

7  X  10-1^ 

1.10 

57.0 

oo 

X 

o 

! 

1.26 

53.1 

9  X  10-^'^ 

1.42 

51.5 

1  X  10-1® 

1.58 

50.1 

ratios  which  are  a  factor  of  1.5  to  4  larger  than  those  provided  by  conventional  Hartmann  sensor  processing. 

The  radially  averaged  intensity  patterns  shown  in  Fig.  10  were  computed  by  averaging  the  beacon  plane 
intensity  for  the  projected  beam,  and  then  averaging  this  result  over  circles  of  constant  radius.  These  plots 
indicate  that  the  nonlinear  optimization  approach  provided  results  superior  to  the  conventional  Hartmann 
sensor-based  approach  for  all  cases  of  saturated  scintillation  in  both  on-axis  intensity,  and  by  providing 
narrower  beam  width.  In  fact,  the  nonlinear  optimization  technique  for  wave  front  control  yields  results 
which  approach  those  provided  by  the  idealized  branch  point  reconstruction  technique. 
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Strehl  Ratio 


Figure  9:  Strehl  ratio  results  for  the  various  correction  schemes.  In  the  legend  “BP”  refers  results  obtained 
for  (j)Qp  (i,  j),  “LS”  refers  to  results  obtained  for  j),  “WFS”  refers  to  results  obtained  with 

and  “NL  OPT”  refers  to  results  obtained  with  (p^¥{i,j). 


RaiU  Avaragt  Inisnsites 


Figure  10:  Example  radially  averaged  beacon  plane  intensities  for  the  projected  beam:  (a)  =  2  x  10“^^, 

<T^  =  0.316;  (b)  =  4  X  10“^'^,  =  0.631;  and  (c)  =  7  x  10“^'^,  =  1.104.  In  the  legend  “BP” 

refers  results  obtained  for  4>gp  {i,j),  “LS”  refers  to  results  obtained  for  <pL^{i,j),  “WFS”  refers  to  results 
obtained  with  and  “NL  OPT”  refers  to  results  obtained  with 
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5  Experiment  and  experimental  results. 


An  experiment  has  been  set  up  and  designed  to  verify  the  ability  to  control  the  amplitude  of  the  transmitted 
wave  front  field  falling  on  the  DM2  in  the  two-deformable-mirror  concept.  In  the  experiment,  a  collimated 
laser  beam  illuminates  a  deformable  mirror  which  was  controlled  electronically.  A  wave  front  sensor  (WFS) 
measurement  of  the  wave  front  field  reflected  from  the  deformable  mirror  is  taken  by  a  digital  camera.  The 
WFS  measurement  is  used  to  obtain  the  phase  of  the  reflected  wave  front  field.  A  Fourier-transforming 
lens  is  used  to  obtain  the  Fraunhofer  diffraction  or  far  field  intensity  of  the  reflected  wave  front  field.  The 
Fraunhofer  diflFraction  pattern  of  the  reflected  wave  front  field  is  the  amplitude  of  the  transmitted  wave 
front  field  falling  on  DM2  in  the  two-deformable-mirror  concept.  The  far  field  intensity  of  the  reflected 
wave  front  field  is  measured  by  a  digital  camera.  The  far  field  intensity  pattern  is  compared  with  the 
theoretical  far  field  pattern.  The  theoretical  far  field  intensity  pattern  is  obtained  by  taking  the  Fourier 
transform  of  exp{jipi{x)},  where  rpi{x)  represents  the  reconstructed  phase  from  the  WFS  measurement. 
In  the  experiment,  a  good  match  was  obtained  between  the  measured  intensity  pattern  and  the  intensity 
pattern  predicted  from  the  Hartmann  sensor  measurement. 

The  nonlinear  response  and  strong  coupling  of  control  channels  in  the  deformable  mirror  make  the 
deformable  mirror  used  here  difficult  to  drive  to  a  desired  surface  shape  [33,  34].  Each  actuator  produces 
an  electrostatic  force  to  push  or  pull  the  mirror  surface  to  form  a  particular  mirror  surface  shape,  and  each 
actuator  modulates  the  whole  surface  of  the  mirror.  Therefore,  no  linear  relationship  between  the  control 
signals  and  the  exact  deformation  of  the  mirror  exists. 

The  deformable  mirror  used  in  the  experiment  is  the  micromachined  membrane  deformable  mirror 
(MMDM)  from  OKO  technologies  [34,  35].  MMDM  consists  of  silicon  chip  mounted  over  a  PCB  holder. 
The  chip  contains  silicon  nitride  membrane,  which  is  coated  with  aluminum  to  form  the  mirror.  The  PCB 
contains  the  control  electrode  structure,  spacer,  and  connector.  The  design  of  the  assembled  mirror  is  shown 
in  Fig  11.  The  MMDM  consists  of  37  actuators  which  drive  the  deformable  mirror  surface.  A  layout  of  the 
MMDM  actuators  is  shown  in  Fig.  12.  The  shape  of  the  reflective  membrane  is  controlled  by  applying  volt- 
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ages  between  the  membrane  and  the  control  electrodes.  These  control  voltages  are  supplied  by  8  bit  digital 
boards  installed  in  a  computer  and  two  identical  high-voltage  DC  amplifier  boards.  The  output  voltages 
of  the  digital  board  are  in  the  range  from  1  to  3.65  V.  The  range  of  the  output  voltage  has  been  digitized 
between  0  and  255.  Then,  the  output  voltages  are  sent  to  two  identical  high-voltage  amplifier  boards.  Each 
board  contains  20  non-inverting  DC  amplifiers.  Each  amplifier  has  a  gain  59,  providing  the  output  voltage 
swing  of  1  to  213  V.  The  amplified  output  voltages  are  sent  to  MMDM  and  drive  the  deformable  mirror 
surface. 

The  nonlinear  response  and  strong  coupling  of  control  channels  in  the  deformable  mirror  make  this 
deformable  mirror  difficult  to  drive  to  a  desired  surface  shape  [33,  34].  Each  actuator  will  produce  an 
electrostatic  force  to  push  or  pull  the  mirror  surface.  The  force  produce  by  each  actuator  is  defined  as: 

(41) 

where  eco  is  the  dielectric  constant  of  the  material  between  the  electrodes  (air),  5  is  the  area  of  the  electrode, 
d  is  the  distance  between  the  electrodes,  Vt  is  the  bias  voltage  and  Vc  is  the  control  voltage.  Each  actuator 
modulates  the  whole  surface  of  the  mirror.  The  overall  mirror  surface  deformation  has  a  nonlinear  response. 
The  strong  coupling  of  control  channel  is  presented  because  the  actuators  are  close  to  each  other,  and  due 
to  the  mechanical  properties  of  the  surface  material.  Therefore,  we  were  unable  to  estimate  the  relationship 
between  the  control  signals  and  the  exact  deformation  of  the  mirror.  The  control  signals  can  be  found  from 
the  wave  front  phase  by  the  least  squares  phase  reconstruction  technique  [3]. 

An  influence  function  and  the  corresponding  WFS  image  are  shown  in  Fig.  13.  The  influence  function 
is  obtained  by  applying  voltage  level  255  to  actuator  2  and  64  to  other  actuators.  During  the  phase  recon¬ 
struction,  only  center  980x980  pixels  are  used  from  the  WFS  image.  Also,  the  influence  function  1,  19  and 
37  are  shown  in  Fig.  14.  A  calibration  phase  is  obtained  by  replacing  the  deformable  mirror  with  the  fiat 
mirror.  This  calibration  phase  is  the  deterministic  phase  error  caused  by  the  deterministic  tilt  in  the  wave 
front  field,  and  is  subtracted  from  each  reconstructed  phase.  The  calibration  phase  and  the  corresponding 
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Table  2:  Focal  lengths  of  lenses  used  in  the  experiment  setup. 


Lens  label 

Focal  length  (mm) 

Lo 

400 

Li 

175 

Lft 

500 

Ltl 

75.6 

Table  3:  Nominal  distances  used  in  the  experiment  setup. 


Element  Separation 

Iris  to  Center  of  Beam  Splitter  I 

145 

Center  of  Beam  Splitter  I  to  Center  of  Beam  Splitter  II 

Center  of  Beam  Splitter  II  to  Deformable  Mirror 

Center  of  Beam  Splitter  II  to  Lo 

173 

Lq  to  Li 

Li  to  Hartmann  Lenslet  Array 

193 

Hartmann  Lenslet  Array  to  Hamamatsu  Digital  Camera 

26 

Center  of  Beam  Splitter  I  to  Lft 

97 

Lft  to  Turning  Mirror 

330 

Turning  Mirror  to  Ltl 

Ltl  to  Cohu  Digital  Camera 

WFS  image  are  shown  in  Fig.  15. 

In  order  to  apply  the  phase  deformation  to  DMi  and  DM2,  a  set  of  actuator  commands  are  sent  to  the 
DMi  and  DM2.  The  actuator  commands  can  be  obtained  by  using  the  least  squares  phase  reconstruction 
technique  [3,  pages  177-197].  The  influence  function  is  obtained  by  applying  the  voltage  level  255  to  the 

actuator  and  64  to  other  actuators. 

An  experiment  was  designed  to  verify  the  ability  to  control  the  amplitude  of  the  wave  front  field  falling  on 
the  DM2  in  the  two-deformable-mirror  concept.  The  block  diagram  of  the  experimental  setup  for  controlling 
the  amplitude  falling  on  the  DM2  is  shown  in  Fig.  16.  Lens  focal  lengths  and  nominal  positioning  dimensions 
in  the  experiment  setup  are  summarized  in  Table  2  and  3. 

A  17  mW  He-Ne  laser  beam  with  A=632.8  nm  was  used  as  a  light  source.  The  laser  beam  passed  through 
the  neutral-density  filter.  The  neutral-density  filter  is  used  to  reduce  the  intensity  of  the  laser  beam.  The 
laser  beam  after  the  neutral-density  filter  passed  through  a  beam  expander.  The  purpose  of  the  beam 
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expander  is  to  generate  an  expanded  and  collimated  laser  beam.  Then,  the  expanded  and  collimated  laser 
beam  passed  through  an  iris  to  reduce  the  diameter  of  the  laser  beam  to  15  mm.  The  laser  beam  after 
the  iris  illuminated  on  the  Beam  Splitter  I.  The  Beam  Splitter  I  directed  the  laser  beam  to  Beam  Splitter 
II.  The  Beam  Splitter  II  directed  the  leiser  beam  to  the  deformable  mirror.  A  wave  front  field  is  reflected 
from  the  deformable  mirror.  The  deformable  mirror  is  the  micromachined  membrane  deformable  mirror 
(MMDM)  from  OKO  Technologies  [35].  The  reflected  wave  front  field  is  reflected  back  to  Beam  Splitter 
11.  Beam  Splitter  II  sends  half  of  the  wave  front  toward  the  afocal  telescope  and  half  of  the  wave  front  to 
Beam  Splitter  I.  The  afocal  telescope  is  formed  by  two  convex  lens.  The  purpose  of  the  afocal  telescope 
is  to  reduce  the  diameter  of  the  beam  from  15mm  to  6.56mm.  The  mathematical  expression  for  the  beam 
diameter  reduction  caused  by  an  afocal  telescope  is  represented  as  follow: 


where  ^2  is  the  diameter  of  the  wave  front  field  after  the  afocal  telescope,  di  is  the  diameter  of  the  wave  front 
field  before  the  afocal  telescope,  /i  is  the  focal  length  of  the  lens  Lq  and  /2  is  the  focal  length  of  the  lens  Li. 
In  the  experiment,  d2  is  6.56  mm,  di  is  15  mm,  /i  is  400  mm  and  /2  is  175  mm.  The  compacted  wave  front 
field  is  passed  to  the  lenslet  array.  The  lenslet  array  segments  the  wave  front  field  and  focusses  the  light  on 
the  detector  plane  of  the  Hamamatsu  digital  camera.  Each  lenslet  of  the  lenslet  array  has  focal  length  =  24 
mm  and  side  length  =  328  ^m.  Hence,  there  were  20  subapertures  across  the  beam  diameter.  The  model 
number  of  the  Hamamatsu  digital  camera  is  C4742-95.  The  detector  plane  of  the  Hamamatsu  digital  camera 
consists  of  a  1024x1024  pixel  array,  and  the  rectangular  pixels  have  dimensions  of  6.7  /imx6.7  fim.  A  WFS 
image  is  taken  by  the  Hamamatsu  digital  camera.  The  WFS  image  is  processed  by  the  phase  reconstruction 
software  to  reconstruct  the  wave  front  phase.  The  theoretical  far  field  intensity  can  be  computed  by  the 
following  expression 

ItH  =  |FT[exp{#i(x)}]l2  (43) 

where  Fr[']  is  the  Fourier  transform  operator  and  tpiix)  is  the  reconstructed  phase. 
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The  dynamic  range  of  Fourier  spectrum  predicted  by  Eq.  (43)  is  higher  than  the  dynamic  range  of  the 
detector  in  the  Cohu  digital  camera  [36,  pages  92-93].  Hence,  only  the  brightest  part  of  the  far  field  intensity 
which  computed  from  Eq.  (43)  is  visible.  In  order  to  see  the  far  field  intensity  in  detail,  the  theoretical  far 
field  intensity  can  be  computed  by  the  mathematical  function  [36], 

/e,  =  log(/,„  -h  10000)  (44) 


The  logarithm  function  is  used  to  reduce  the  dynamic  range  of  the  far  field  intensity. 

The  other  half  of  the  wave  front  field  from  Beam  Splitter  II  passed  through  Beam  Splitter  I  and  fell 
on  the  lens  LpT  with  /=500  mm.  This  lens  is  the  Fourier-transforming  lens  in  the  two-deformable-mirror 
concept.  The  wave  front  field  after  the  lens  was  passed  toward  the  turning  mirror.  The  reason  for  using  the 
turning  mirror  is  to  create  more  space  on  the  optical  bench.  The  wave  front  field  reflected  from  the  turning 
mirror  is  propagated  through  two  neutral- density  filters  to  avoid  saturation  the  detector.  A  positive  lens 
Ltl  with  /=75.6  mm  is  used  to  magnify  the  far  field  image.  The  ratio  of  the  magnification  is  3:1.  The  thin 
lens  law  is  used  to  achieve  the  magnification.  Prom  the  thin  lens  law, 


M  =  - 


£i 

) 

So 


(45) 


and 


(46) 


where  sign  indicate  the  image  is  inverted,  M  is  the  magnification,  Si  is  the  image  distance,  Sq  is  the  object 
distance  and  /  is  the  focal  length  of  the  lens  which  used  in  the  magnification  system.  In  the  experiment,  / 
is  75.6mm.  The  Si  and  So  cannot  be  measured  exactly  from  the  experiment  because  the  plane  of  the  lens 
cannot  be  located.  A  variable  neutral-density  filter  wheel  is  placed  in  front  of  the  Cohu  digital  camera  to 
control  exposure  on  the  detector  plane  of  the  Cohu  digital  camera,  which  had  a  fixed  integration  time  of  33 
ms.  The  model  number  of  the  Cohu  digital  camera  is  4110.  The  image  plane  has  480x752  pixels,  and  the 
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rectangular  pixels  have  dimensions  of  8.66  /imx9.92  /im. 

The  location  of  focus  of  both  digital  cameras  can  be  found  by  counting  the  number  of  pixels  between  the 
first  zeros  of  the  point  spread  functions  with  the  deformable  mirror  replaced  by  a  flat  mirror.  The  deformable 
mirror  is  not  used  for  alignment  because  deformable  mirror  surface  is  not  flat.  The  deformable  mirror  has 
initial  deviation  from  plane  of  less  than  0.6  /xm  [35].  In  the  case  of  finding  the  location  of  focus  of  Cohu 
digital  camera,  the  diameter  of  the  laser  beam  was  reduced  to  8mm  from  15  mm  to  broaden  the  Airy  disk, 
making  it  easier  to  find  focus.  The  width  between  the  first  zeros  of  the  point-spread  functions  on  the  Cohu 
digital  camera  can  be  calculated  by  the  following  expression: 

2.44M/A 

width  = - (47) 

where  M  is  the  magnification  of  the  image,  /  is  the  focal  length  of  the  Fourier-transforming  mirror,  A  is  the 
wavelength  of  the  laser  beam  and  D  is  the  diameter  of  the  laser  beam.  In  the  experiment,  M=3,  /=500  mm, 
A=632.8  nm  and  D=8  mm.  Therefore,  the  width  between  the  first  zeros  of  the  point-spread  functions  on 
the  Cohu  digital  camera  is  289.51  /xm.  The  number  of  pixels  which  span  the  x-axis  slice  of  the  point  spread 
function  on  the  detector  plane  is  thus  33  pixels.  From  the  experiment,  we  found  that  the  number  of  pixels  in 
the  x-axis  slice  of  the  point  spread  function  is  34  pixels.  The  picture  and  the  x-axis  slice  of  the  point  spread 
function  are  shown  in  Fig.  17  and  Fig.  18.  Finally,  the  far  field  intensity  is  measured  and  displayed  on  the 
Cohu  digital  camera. 

In  order  to  test  the  performance  of  the  experiment,  different  sets  of  micromachined  membrane  deformable 
mirror  actuator  voltage  signals  are  applied  to  drive  the  mirror  to  different  surface  shapes.  In  a  particular  set 
of  mirror  actuator  voltage  signals,  a  WFS  image  on  the  Hamamatsu  digital  camera  and  a  far  field  intensity 
on  the  Cohu  digital  camera  are  measured.  The  phase  of  the  wave  front  field  reflected  from  the  deformable 
mirror  is  reconstructed  from  the  wave  front  sensor  measurement.  The  measured  far  field  intensity  pattern 
can  be  compared  with  the  theoretical  intensity  far  field  pattern.  In  order  to  compare  the  measured  far  field 
intensity  pattern  and  the  theoretical  far  field  intensity  pattern,  same  spacing  in  both  image  are  required. 
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The  measured  far  field  intensity  pattern  has  been  interpolated  by  using  the  bicubic  interpolation  to  have 
same  spacing  with  the  theoretical  far  field  intensity  pattern.  Both  far  field  intensity  patterns  have  256x256 
pixel  array,  and  the  rectangular  pixels  have  dimensions  of  48.04  //mx 48.04  /^m. 

The  experiment  results  shown  in  Fig.  20  are  obtained  by  applying  the  maximum  voltage  (in  digitized 
levels)  of  255  to  actuators  one  and  two,  and  constant  voltage  64  on  the  other  actuators.  Figure  20(a)  show 
the  reconstructed  phase.  Fig.  20(b)  shows  the  theoretical  far  field  intensity  pattern,  and  Fig.  20(c)  shows 
the  measured  far  field  intensity  pattern  on  the  Cohu  digital  camera. 

The  experiment  results  that  shown  in  Fig.  21  is  obtained  by  applying  the  maximum  voltage  of  255  (in 
digitized  levels)  to  actuator  1,  19  and  37,  and  the  constant  voltage  64  to  the  other  actuators.  Fig.  21(a) 
shows  the  reconstructed  phase,  Fig.  21(b)  shows  the  theoretical  far  field  intensity  pattern,  and  Fig.  21(c) 
shows  the  measured  far  field  intensity  pattern  on  the  Cohu  digital  camera. 

From  the  experiment  results  above,  the  far  field  intensity  pattern  is  successfully  obtained.  This  far  field 
intensity  is  the  intensity  of  the  transmitted  wave  front  field  falling  on  the  DM2  in  the  two-deformable-mirror 
concept.  The  ringing  effect  present  in  all  theoretical  far  field  intensity  patterns  is  due  to  the  two-dimensional 
triangle  function  which  were  used  as  elementary  functions  in  phase  reconstruction. 

6  Conclusions 

In  this  program  the  problem  of  correcting  for  atmospheric  turbulence-induced  scintillation  errors  in  laser 
beam  projection  and  imaging  through  the  atmosphere  has  been  addressed.  The  approach  taken  was  the 
implementation  and  study  of  a  phase  and  amplitude  correction  technique  which  uses  two  deformable  mir¬ 
rors,  and  then  extend  these  results  to  look  at  a  novel  algorithm  for  controlling  a  single  deformable  mirror. 
This  work  has  been  successful  in  that  it  has  been  demonstrated  that  the  possibility  of  compensating  for 
turbulence-induced  scintillation  exists  for  both  the  laser  beam  projection  and  imaging  applications.  Studies 
of  the  implications  of  using  real  adaptive  optical  system  hardware  to  implement  the  two  deformable  mirror 
technique  have  shown  the  limitations  of  this  approach.  However,  the  results  obtained  here  for  the  nonlinear 
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optimization-based  approach  to  controlling  a  single  deformable  mirror  show  a  great  deal  of  promise,  and  will 
be  pursued  in  the  future. 
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Figure  11:  The  schematic  diagram  of  the  micromachined  membrane  deformable  mirror. 
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WFS  Image 
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Figure  13:  (a)  WFS  image  (b)  The  influence  function;  with  applying  voltage  level  255  to  actuator  2  and  64 
to  other  actuators  of  the  deformable  mirror. 
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tnfluanca  function  Hh 


Figure  14;  The  influence  functions  are  obtained  by  applying,  maximum  voltage  255  on  actuator  and 
constant  voltage  64  on  the  other  actuators:  (a)  The  influence  function  1.  (b)  The  influence  function  19.  (c) 
The  influence  function  37. 
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WFS  Image  with  the  flat  mirror 
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Figure  15;  (a)  WFS  image  (b)  The  calibration  phase;  with  the  deformable  mirror  replaced  by  flat  mirror. 
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Cohu  Digital  Camera 


Neutral  Density 
Filter  Wheel 

Neutral  Density 
Filters 


Turning  Mirror 


Beam  Expander 


Neutral  Density 
Filter 


t 


Laser 


|||  Hamamatsu  Digital  Camera 
®®®@  Hartmann  Lenslet  Array 

A- 

Deformable  Mirror 

Beam  Splitter  II 


Figure  16:  Block  diagram  of  the  experimental  setup  for  controlling  the  amplitude  falling  on  DM2. 
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Figure  17;  The  point  spread  function  with  the  flat  mirror  falling  on  the  Cohu  digital  camera. 
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Figure  19:  Photograph  of  the  experiment  setup. 
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M««su«d  lar  fl«ld  Intensity  pattern 


Figure  20:  The  experiment  results  is  obtained  by  applying  the  maximum  voltage  of  255  (in  digitized  levels) 
on  actuators  one  and  two,  and  the  constant  voltage  of  64  on  the  other  actuators:  (a)  Reconstructed  phase; 
(b)  theoretical  far  field  intensity  pattern;  (c)  measured  far  field  intensity  pattern 


Figure  21:  The  experimental  results  obtained  by  applying,  the  maximum  voltage  of  255  (in  digitized  levels) 
to  actuator  1,  19  and  37,  and  constant  voltage  64  to  the  other  actuators:  (a)  Reconstructed  phase;  (b) 
theoretical  far  field  intensity  pattern;  (c)  measured  far  field  intensity  pattern 
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